Mapas con ggplot

1 Introducción

1.1 Acerca de este documento.

Este documento presenta algunos ejercicios básicos de cartografía digital con R y ggplot2::. En modo representa las capacidades completas de este software para el manejo de sistemas de información geográfica y tampoco se propone presentar resultados pertinentes de investigación. Simplemente presenta algunos objetivos y los pasos técnicos necesarios para llevarlos a cabo.

1.1.1 ¿Por qué ggplot2:: para producir mapas?

Si estamos usando ggplot2:: para el resto de los gráficos de nuestro trabajo mantenernos en el misma librería de funciones para los mapas tiene mucho sentido.

• Nos permite mantener un estilo visual unificado con poco esfuerzo: todos los gráficos -incluyendo mapas- tendrán la misma tipografía, tamaño relativo, convención de títulos y viñetas, colores o formas, etc.

• Nos permite sacar provecho de una sintaxis que ya conocemos. Al producir objetos de la clase ggplot2 podemos modificar muchos aspectos de estos mapas con la misma sintaxis que usamos para cualquier otro gráficos. Las funciones labs() y annotate()producen anotaciones, scale_fill_*() nos permite usar paletas personalizadas de colores, etc.

• Nos permite aprovechar algunas extensiones generales de R y aplicarlas a mapas. Por ejemplo, ggreppel() para evitar la superposición de etiquetas.

## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats
## lag():    dplyr, stats

Temas cubiertos en este documento.

Tema	Métodos	Librería y función
Importación de datos	Desde shapefiles .shp	rgdal::readOGR()
	Raster desde la API de GoogleMaps	ggmap::get_map()
	Rutas desde un teléfono celular	plotKML::readGPX()
Mapas coropléticos	Con polígonos primitivos	ggplot2:: geom_polygon()
	Con mapas+datos	ggplot() + geom_map()
	Con una función ad hoc	mxmaps::
Etiquetado de mapas	Como texto	ggplot2::geom_label
	Ubicación en centroides	sp::coordinates()
	Etiquetado parcial	ifelse() + ggreppel::geom_label_repel()

1.2 Nota sobre la información geográfica del INEGI y los sistemas de coordenadas.

El INEGI tiene puestos a disposición del público shapefiles con distintos niveles de información. Sin considerar los históricos esta información se regular por el Marco Geoestadístico Censal, cuya última versión es de 2016. Incluye las siguientes divisiones y subdivisiones del territorio:

Subdivisiones territoriales básicas del INEGI.

Código	División	Unidades territoriales
AGEE	Estatales	32
AGEM	Municipales	2458
AGEB	Censal	17470
PLUR	Localidades	49498

Nota aparte merece el sistema de coordenadas que utiliza INEGI, también definido en el el Marco Geoestadístico Censal 2016. La definición técnica es MEXICO_ITRF_2008_LCC, con proyección Cónica Conforme de Lambert (CCL) y datum es ITRF2008. Se trata de un sistema de coordenadas rectangular, cuya unidad de medida es el metro. Esta cuadricula cubre el espacio entre los paralelos 17.5N y 29.5N con el meridiano 102 como centro. En términos prácticos esto tiene algunas implicaciones:

• La cartografía que generemos a partir de shapefiles del INEGI se expresa directamente en una proyección ad hoc para México. Por lo tanto:
  • Podemos esperar mapas con pocas deformaciones.
  • Al utilizar unidades rectas para expresar la distancia (y no grados, minutos, segundos, unidades de arco) es fácil interpretar la escala de nuestros gráficos. Esta es razonablemente homogénea.
    
  • La mayor parte de la información geográfica que producen las instituciones gubernamentales de México (INEGI, CONAPO, CONABIO, etc.) utiliza este mismo marco de referencia, por lo que podemos combinar estas capas sin mayores inconvenientes. El INE es un caso aparte, su sistema de información geográfica es diferente.
• Sin embargo otras capas de información geográfica pudiéramos aprovechar y que se produce fuera de la órbita de estas instituciones suele estar expresada en otro sistema de coordenadas. Por ejemplo, en unidades hexadecimales de arco, es decir, latitud y longitud medidas en grados, minutos y segundos. Por lo tanto:
  • No podemos unir esas capas sin una transformación previa, de lo contrario la ubicación de los polígonos no va a ser compatible.
    • Debemos llevar la información no INEGI al sistema de INEGI o viceversa.
      • INEGI dispone de una app en línea para convertir información en lat- long a lcc en http://www.inegi.org.mx/geo/contenidos/geodesia/traninv.aspx. Hasta ahora no he visto una opción para convertir una base de puntos datos completa y mucho menos polígonos. Debemos transformar punto por punto.
      • No encontré una función o librería en R que nos permita hacer este traspaso de manera simple, aunque no he sido exhaustivo en la búsqueda.
    • http://www.gadm.org/ dispone de shapefiles y objetos geoespaciales SPDF de R (.rds) de México, subdivididos hasta nivel municipal con un sistema de coordenadas geográficas. Se trata de archivos pequeños que se cargan directamente en R sin necesidad de importación. Sin embargo no dispone de unidades territoriales menores al municipio, red caminera, áreas y recursos naturales ,etc. y no es fuente oficial. Tampoco podemos estar completamente seguros de su actualización. Úselo bajo su propia responsabilidad.
      
  • La cartografía de GoogleMaps está disponible a través de una API y R tiene varios clientes para esa API. El problema señalado de los diferentes sistemas de coordenadas dificulta su aprovechamiento.
  • R tiene la librería mapproj:: para proyección de mapas y ggplot nos permite utilizarla fácilmente con coord_map(). Sin embargo esta librería sólo admite coordenadas geográficas, por lo que no podemos utilizarla con datos crudos del INEGI.

2 Mapas coropléticos de México con R.

Los mapas coropléticos representan polígonos de unidades territoriales a los que se asigna un color de acuerdo al valor de una variable pertinente. Este color puede representar tanto una escala continua -en la que se genera un degradado continuo de color- o una escala discreta -en la que cada categoría de esa variable tiene un color asignado.

El procedimiento en R usando ggplot2:: es el siguiente:

1. Importar los polígonos con los contornos de las unidades territoriales. Esto producirá un objeto de la clase SPDF o Spatial Polygon Data Frame.

• Los objetos SPDF no son directamente un data.frame, en sentido estricto pertenecen a la clase list. Sin embargo se transforman fácilmente a data.frame.
• Fueron incorporados en R a partir de la librería sp::, un conjunto de clases y métodos para datos espaciales. sp:: es la librería primitiva para información geográfica en R, prácticamente todas las demás depende de un modo u otro de las especificaciones, funciones y métodos de sp::.

2. Unir los datos estadísticos que controlarán el color con la estructura de datos poligonal.

• A través de un join o con el método simplificado de geom_map().

3. Generar el mapa con una llamada a ggplot() y los elementos geométricos geom_polygon() o geom_map().

2.1 Unidad territorial: Estados.

Vamos a utilizar al Estado como unidad territorial, para mantener las cosas simples y reducir el tiempo de procesamiento.

-Instalamos la librería rgdal, que tiene la función readOGR para cargar shapefiles en R y convertirlos en un objeto espacial de R.

• Descargamos desde acá http://www.conapo.gob.mx/es/CONAPO/Datos_Abiertos_del_Indice_de_Marginacion los datos de marginación del CONAPO a nivel municipal. Se pueden agregar por Estados.

• Descargamos los shapefiles del Marco Geoestadístico 2016 del sitio del INEGI (202MB) desde este enlace: http://internet.contenidos.inegi.org.mx/contenidos/productos//prod_serv/contenidos/espanol/bvinegi/productos/geografia/marc_geo/702825217341_s.zip

• Descomprimimos todos los archivos del .zip en el directorio de trabajo.
  • Usamos getwd() para saber cuál es nuestro directorio de trabajo.

2.1.1 Importación de capas y mapas simples.

library(rgdal)

## Loading required package: sp

## rgdal: version: 1.2-5, (SVN revision 648)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/mpaladino/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/mpaladino/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-4

#Cargo el shapefile de los estados y lo asigno a capa_estados.
#"." es el directorio, layer="Estados" el archivo .shp, pero sin la extensión, ya que readOG

capa_estados <-readOGR(".", layer="Estados")

## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "Estados"
## with 32 features
## It has 5 fields

#¿Qué habrá acá dentro?

class(capa_estados) 

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

#str(capa_estados) Demasiado largo, supera a la consola. 

capa_estados$NOM_ENT #Con sensacionales errores de codificación!

##  [1] Distrito Federal                Guerrero                       
##  [3] México                         Morelos                        
##  [5] Sinaloa                         Baja California                
##  [7] Sonora                          Baja California Sur            
##  [9] Zacatecas                       Durango                        
## [11] Chihuahua                       Colima                         
## [13] Nayarit                         Michoacán de Ocampo           
## [15] Jalisco                         Chiapas                        
## [17] Tabasco                         Oaxaca                         
## [19] Guanajuato                      Aguascalientes                 
## [21] Querétaro                      San Luis Potosí               
## [23] Tlaxcala                        Puebla                         
## [25] Hidalgo                         Veracruz de Ignacio de la Llave
## [27] Nuevo León                     Coahuila de Zaragoza           
## [29] Tamaulipas                      Yucatán                       
## [31] Campeche                        Quintana Roo                   
## 32 Levels: Aguascalientes Baja California Baja California Sur ... Zacatecas

capa_estados$CVE_ENT #Este es el que sirve. 

##  [1] 09 12 15 17 25 02 26 03 32 10 08 06 18 16 14 07 27 20 11 01 22 24 29
## [24] 21 13 30 19 05 28 31 04 23
## 32 Levels: 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 ... 32

#A ver...
ggplot() +  
  geom_polygon(data=capa_estados, aes(x=long, y=lat, group=group))

## Regions defined for each Polygons

#Latitudd y longitud Marco Geográfico INEGI.

#Con los contornos de los estados. 
ggplot() +  
  geom_polygon(data=capa_estados, aes(x=long, y=lat, group=group), 
               fill="white", color="black") #Simplemente pasé las líneas a negro y el relleno a blanco fuera de aes(). 

## Regions defined for each Polygons

Obtuvimos un mapa completo de México. Valdría la pena hacer un mapa coroplético en el que coloreamos a cada Estado por la media de Marginación Municipal.

2.1.2 Colorear los Estados por la media del IM. Con join. Más complejidad, más control.

El elemento geométrico primitivo que ggplot2 utiliza para generar mapas a partir de polígonos –como los que importamos desde shapefiles– es geom_polygon(). El data.frame resultante de la conversión de un objeto SPDF es sumamente largo, tiene una fila por cada cambio en un polígono y un columna de grupo para separar cada polígono. Para un mapa coroplético podríamos usar el argumento fill= para los polígonos. Antes debemos resolver el problema de las estructuras de datos. Por el momento la información geográfica está en una estructura de lista -no data.frame- y tiene una entrada por cada curva en los polígonos. Es decir, debemos pasar el objeto SPDF a data.frame y encontrar la manera de compatibilizar dos estructuras de datos con largos diferentes: la del mapa con un fila por cada curva de polígono y la de las medias del Índice de Marginación, con 32 filas, una por cada Entidad Federativa. Haremos lo siguiente:

1. Crear un data frame con el sumario estadístico, en este caso las medias por entidad. Deberá incluir una clave en común el data.frame de polígonos.
2. Convertir el objeto SPDF a data.frame.
3. Aplicar un join para unir ambas estructuras de datos. El dato de interés que controlará el color se repetirá una vez por cada fila, de modo que los largos sean iguales. dplyr::inner_join() se encarga de hacerlo.
4. Producir el mapa con el elemento geométrico geom_polygon(), agregando un argumento fill()

#Cargo la base de marginación y la limpio. 
marginacion <- read_csv("C:/Users/mpaladino/Dropbox/martinpaladino.github.io/datos/Base_Indice_de_marginacion_municipal_90-15.csv", 
    col_types = cols(`AÑO` = col_character()), 
    locale = locale(encoding = "UTF-8")) %>%  
      mutate (POB_TOT=gsub(" ","", POB_TOT)) %>%
      mutate(POB_TOT=as.integer(POB_TOT)) %>%
      filter(ENT!="Nacional")

## Warning: 21915 parsing failures.
##  row     col expected actual
## 2404 SPRIM   a double      -
## 2404 OVSDE   a double      -
## 2404 VHAC    a double      -
## 2404 OVPT    a double      -
## 2404 PL<5000 a double      -
## .... ....... ........ ......
## See problems(...) for more details.

#Obtengo las media de IM para 2015. 
marginacion %>%
  filter(AÑO=="2015") %>% 
  group_by(CVE_ENT, ENT) %>% 
  summarise (mediaIM=mean(IM)) %>% 
  mutate(id=CVE_ENT) -> #Renombre CVE_ENT a id para que coincida con el nomnre de columna del data.frame de los polígonos. 
  mediaIM

mediaIM

## Source: local data frame [32 x 4]
## Groups: CVE_ENT [32]
## 
##    CVE_ENT                  ENT    mediaIM    id
##      <chr>                <chr>      <dbl> <chr>
## 1       01       Aguascalientes -0.9220909    01
## 2       02      Baja California -1.5100000    02
## 3       03  Baja California Sur -1.2734000    03
## 4       04             Campeche -0.1338182    04
## 5       05 Coahuila de Zaragoza -1.1691579    05
## 6       06               Colima -1.0235000    06
## 7       07              Chiapas  0.8246441    07
## 8       08            Chihuahua -0.3207761    08
## 9       09     Distrito Federal -1.7696875    09
## 10      10              Durango -0.3135128    10
## # ... with 22 more rows

#Un df de 32x4 con "id" igual al SPDF.

#Paso el objeto SPDF a data.frame con `fortify()`, una función de ggplot2 con métodos para convertir objetos diversos a data.frame.

capa_estados_df <- fortify(capa_estados, region="CVE_ENT") #region="CVE_ENT", el agrupamiento de los polígonos. Coincide en nombre y contenido con mediaIM$id.    

capa_estados_df_mediaIM <- inner_join(capa_estados_df, mediaIM, by="id") #Uno las medias con los polígonos por la variable clave: id.

head(capa_estados_df_mediaIM)   #La media se repite una vez por punto de polígono, unidas por id. 

##      long     lat order  hole piece id group CVE_ENT            ENT
## 1 2470518 1155029     1 FALSE     1 01  01.1      01 Aguascalientes
## 2 2470552 1154983     2 FALSE     1 01  01.1      01 Aguascalientes
## 3 2470607 1154988     3 FALSE     1 01  01.1      01 Aguascalientes
## 4 2470637 1155017     4 FALSE     1 01  01.1      01 Aguascalientes
## 5 2470661 1155046     5 FALSE     1 01  01.1      01 Aguascalientes
## 6 2470683 1155073     6 FALSE     1 01  01.1      01 Aguascalientes
##      mediaIM
## 1 -0.9220909
## 2 -0.9220909
## 3 -0.9220909
## 4 -0.9220909
## 5 -0.9220909
## 6 -0.9220909

ggplot(capa_estados_df_mediaIM) +  
  geom_polygon(aes(x=long, y=lat, 
                   group=group,      #El argumento group=group arma grupos tanto para para los polígonos como para `fill=`. 
                   fill=mediaIM))    #Aquí especifico la columna que controla los colores. 

2.1.3 Colorear los Estados por la media del IM. Con geom_map()

ggplot2 tiene un elemento geométrico específico para generar mapas. Usa el primitivo geom_polygon(), como lo hicimos en el ejercicio anterior, pero permite especificar directamente el mapa y automatiza el manejo de múltiples estructuras de datos. Por lo tanto no es necesario hacer unir previamente las estructuras de datos. Esto simplifica la

#Crear un data.frame con los polígonos, incluyendo "CVE_ENT"
#Retomo las estructuras de datos que ya cree: capa_estados_df y mediaIM. 

#Sintaxis de geom_map.

ggplot(mediaIM,                          #Los datos que pasamos a ggplot son el resumen estadístico de intereés, NO el mapa.  
       aes(map_id = id)) +               #map_id es la columna con la clave de unión. Nos ahorra el join previo. 
    geom_map(aes(fill = mediaIM),        #fill= columna que controla el color.
             map = capa_estados_df) +    #map= objeto SPDF convertido a data.frame. 
  expand_limits(x = capa_estados_df$long, y = capa_estados_df$lat)  #Especificamos el tamaño del mapa al máximo de lat y long.

Obtuvimos el resultado que buscábamos y el un mapa exploratorio aceptable. Sin embargo podemos mejorarlo usando la sintaxis usual de ggplot2. Esta es la gran ventaja de usar una misma librería para todos nuestros gráficos: conociendo una sintaxis podemos modificar muchos tipos de gráficos, incluyendo mapas. Además podemos aplicar el mismo estilo visual. En este caso usaremos dplyr:: y ggplot2:: para ajustar el gráfico.

1. Creando datos más pertinentes. Media de IM municipal ponderada por población.
2. Agregando una escala de distancia. Aprovechamos la proyección definida por INEGI, que en para el territorio mexicano es plana y está expresada en unidad de medida métrica. Esto no aplica de manera directa para cartografía generada a partir de coordenadas geográficas.
   
3. Ubicando la leyenda de colores, con el título bien formateado, en el Golfo de México, para aprovechar mejor el espacio.
4. Usando el código de colores semáforo, donde rojo marca una advertencia.

marginacion %>%
  filter(AÑO=="2015") %>% 
  group_by(CVE_ENT, ENT) %>% 
  summarise (`Media IM \nPonderada por\npoblación` =     #Uso \n para insertar directamente saltos de línea
               weighted.mean(IM, POB_TOT)) %>%           #weighted.mean() Media ponderada
  mutate(id=CVE_ENT) ->
  mediaIM

ggplot(mediaIM, aes(map_id = id)) + 
  geom_map(aes(fill = `Media IM \nPonderada por\npoblación`), map = capa_estados_df) + 
  expand_limits(x = capa_estados_df$long, y = capa_estados_df$lat) +
  scale_fill_continuous(low="green", high="red") +                           #Cambio la escala contínua de colores.
  labs(title="Media del Índice de Marginación Municipal por Estados 2015",   #Las anotaciones con la sintaxis usual de ggplot2::
       subtitle="Ponderado por población", 
       caption="Elaboración propia. \n Datos geográficos: INEGI 2016. \n Datos de estadísticos: CONAPO 2015") +
  theme_void() +                                                             #Elimina escalas, marcas de coordenadas, etc. 
  geom_errorbarh(aes(x=1.5e6, xmin=1.5e6-1e5, xmax=1.5e6+1e5, y=1e6), height=5e4) + #Uso una barra de error horizontal para la escala
  annotate("text", x= 1.5e6, y=1.0e6+8e4, label="200km", size=2) +           #Nota de la escala +- 100000 metros=200km. 
  theme(legend.position = c(0.8, 0.7))                                       #Ubico la leyenda dentro del gráfico en una zona vacia.

2.2 Con coordenadas geográficas de http://www.gadm.org/.

Global Administrative Areas tiene objetos SPDF –formato nativo de R– con polígonos administrativos de prácticamente todos los países del mundo con sistema de coordenadas geográficas . Al estar en formato de R no es necesario importarlo, simplemente los cargamos a R. El proceso es más rápido.

#Importo y preparo los datos
#===========================
library(mapproj)

## Loading required package: maps

## 
## Attaching package: 'maps'

## The following object is masked from 'package:purrr':
## 
##     map

estados_poly <- readRDS("MEX_adm1.rds") #Directo a R. 
estados_df <- fortify(estados_poly)     #Lo paso a data.frame

## Regions defined for each Polygons

marginacion <- read_csv("C:/Users/mpaladino/Dropbox/martinpaladino.github.io/datos/Base_Indice_de_marginacion_municipal_90-15.csv", 
    col_types = cols(`AÑO` = col_character()), 
    locale = locale(encoding = "UTF-8")) %>%  
      mutate (POB_TOT=gsub(" ","", POB_TOT)) %>%
      mutate(POB_TOT=as.integer(POB_TOT)) %>%
      filter(ENT!="Nacional")

## Warning: 21915 parsing failures.
##  row     col expected actual
## 2404 SPRIM   a double      -
## 2404 OVSDE   a double      -
## 2404 VHAC    a double      -
## 2404 OVPT    a double      -
## 2404 PL<5000 a double      -
## .... ....... ........ ......
## See problems(...) for more details.

marginacion %>%
  filter(AÑO=="2015") %>% 
  group_by(CVE_ENT, ENT) %>% 
  summarise (mediaIM=mean(IM)) %>% 
  mutate(id=as.numeric(CVE_ENT)) -> #Nótese que as.numeric. else los 10 primeros estados no plotean. 
  mediaIM

ggplot(mediaIM, aes(map_id = id)) + 
    geom_map(aes(fill = mediaIM), 
             map = estados_df) + 
    expand_limits(x = estados_df$long, y = estados_df$lat) + 
  theme_minimal() +
  labs(title="Proyección Mercator") #Funciona. Van en serie. 

ggplot(mediaIM, aes(map_id = id)) +   
    geom_map(aes(fill = mediaIM), 
             map = estados_df) + 
    expand_limits(x = estados_df$long, y = estados_df$lat) + 
  coord_map ("cylindrical") + 
  theme_void() + 
  labs(title="Proyeccción cilíndrica")

ggplot(mediaIM, aes(map_id = id)) +   
    geom_map(aes(fill = mediaIM), 
             map = estados_df) + 
    expand_limits(x = estados_df$long, y = estados_df$lat) + 
  coord_map ("gilbert") + 
  theme_void() + 
  labs(title="Proyeccción Gilbert")

ggplot(mediaIM, aes(map_id = id)) +   
    geom_map(aes(fill = mediaIM), 
             map = estados_df) + 
    expand_limits(x = estados_df$long, y = estados_df$lat) + 
  coord_map ("lagrange") + 
  theme_void() + 
  labs(title="Proyeccción Lagrange")

#Con paneles. 

marginacion %>%
  filter(AÑO=="2015") %>% 
  group_by(CVE_ENT, ENT) %>% 
  summarise (`mediaPL>5000`=mean(`PL<5000`), mediaANALF=mean(ANALF), mediaPO2SM=mean(PO2SM), mediaOVSAE=mean(OVSAE)) %>% 
  mutate(id=as.numeric(CVE_ENT)) %>% 
  ungroup () %>% 
  select(-CVE_ENT) %>% 
  gather(., clave, valor, -id, -ENT) ->   mediaIM

ggplot(mediaIM, aes(map_id = id)) +   
    geom_map(aes(fill = valor), map = estados_df) + 
    expand_limits(x = estados_df$long, y = estados_df$lat) + 
  facet_wrap(~clave) +
  scale_fill_continuous(low="green", high="red") + 
  theme_void()

2.3 Unidad territorial: Municipios de Tlaxcala.

#Primer paso: cargar el shapefile de Tlaxcala con polígonos municipales. 

tlaxcala_poly <- readOGR(".", "tlax_municipal")

## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "tlax_municipal"
## with 60 features
## It has 160 fields

tlaxcala_df <- fortify(tlaxcala_poly, region="CVEGEO")

marginacion_tlaxcala <-read_csv("C:/Users/mpaladino/Dropbox/martinpaladino.github.io/datos/Base_Indice_de_marginacion_municipal_90-15.csv", 
    col_types = cols(`AÑO` = col_character()), 
    locale = locale(encoding = "UTF-8")) %>%  
      mutate (POB_TOT=gsub(" ","", POB_TOT)) %>%
      mutate(POB_TOT=as.integer(POB_TOT)) %>%
      filter(AÑO=="2015" & CVE_ENT=="29")

## Warning: 21915 parsing failures.
##  row     col expected actual
## 2404 SPRIM   a double      -
## 2404 OVSDE   a double      -
## 2404 VHAC    a double      -
## 2404 OVPT    a double      -
## 2404 PL<5000 a double      -
## .... ....... ........ ......
## See problems(...) for more details.

names(marginacion_tlaxcala) [1] <- "id" #Coincide con el id de los polígonosPorque filter (o tibble o los dos) complican la cadena con el nombre de la primera variable. Investigar y reportar el bug. si no iría directamente un mutate(), pero falla.

#Validación de polígonos. 
ggplot() +  
  geom_polygon(data=tlaxcala_poly, aes(x=long, y=lat, group=group), 
               fill="white", color="black") + 
  labs(title="Municipios de Tlaxcala", 
       subtitle="Validación de polígonos lat long") +
  theme_minimal() 

## Regions defined for each Polygons

ggplot() +  
  geom_polygon(data=tlaxcala_poly, aes(x=long, y=lat, group=group), 
               fill="white", color="black") + 
  geom_point(aes(
    x=as.data.frame(coordinates(tlaxcala_poly))$V1, #coordinates estima los centroides de cada polígono. Anda mejor en Wyoming...  
    y=as.data.frame(coordinates(tlaxcala_poly))$V2)) + 
  labs(title="Municipios de Tlaxcala", 
       subtitle="Proyección Lagrange") +
  theme_minimal() + coord_map("lagrange")

## Regions defined for each Polygons

#Coroplético por IM. 

ggplot(marginacion_tlaxcala, aes(map_id = id)) +  
  geom_map(aes(fill = IM), 
             map = tlaxcala_df) +  
  expand_limits(x = tlaxcala_df$long, y = tlaxcala_df$lat) + 
  coord_map("lagrange") + 
  theme_minimal() + 
  scale_fill_continuous(low="green", high="red")

#Etiquetado con Nombre de Municipio. 
#La función coordinates() regresa coordenadas x y del centroide de cada polígono, aunque tiene problemas con los polígonos irregulares. Las adosamos a la base con los datos de marginación. 

marginacion_tlaxcala <- cbind(marginacion_tlaxcala, 
                              as.data.frame(coordinates(tlaxcala_poly))) #coordinates() regresa una matríz, tengo que pasarla a df antes de pegarla. Automáticamente les asignará los nombres V1 y V2. 

library(ggrepel)

ggplot(marginacion_tlaxcala, aes(map_id = id)) +  
  geom_map(aes(fill = IM), 
             map = tlaxcala_df) +  
  geom_text_repel(aes(x=V1, y=V2, label=MUN), size=2) +  #Con geom_repel para evitar overplotting. 
  expand_limits(x = tlaxcala_df$long, y = tlaxcala_df$lat) + 
  coord_map("lagrange") + 
  theme_minimal() + 
  scale_fill_continuous(low="green", high="red")

3 Importación de capas adicionales.

Hay muchas formas de importar capas de datos adicionales. Aquí cubriremos una básica: importar directamente de GoogleMaps mapas raster sobre los que podemos graficar polígonos. La función get_map() de ggmap:: permite importar mapas raster de GoogleMaps, OpenStreetMap, Stamen Maps y CloudMade. En general los mejores resultados se obtienen con GoogleMaps, que tiene un API estable y con buenos recursos. OSM tiene problemas frecuentes de no respuesta en los servidores y stamen cambia con frecuencia su API. Debemos especificar las coordenadas que nos interesan como longitud-latitud (ene ese orden), el nivel de zoom (donde 3 es un continente y 21 un edificio), el tipo de mapa (“terrain”, “satellite”, “roadmap”, etc.).

library(ggmap)
google_tlax_terreno <- get_map(
  location=c(-98, 19.4),     #Long y lat del centro del mapa que buscamos
  source="google",           #Fuente, tb OpenStreetView
  maptype="terrain",         #Tipo. También "satellite", "roadmap"
  zoom=9)                    #entre 1 y 21. 1 planisferia, 21 edificio. 

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=19.4,-98&zoom=9&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

#Así lo visualizamos. 
ggmap(google_tlax_terreno)

#Agregamos una capa de divisiones políticas. 
ggmap(google_tlax_terreno) + geom_polygon(data=tlaxcala_poly, aes(x=long, y=lat, group=group), color="black", fill=NA) #Evito que rellene con gris. 

## Regions defined for each Polygons

3.1 Importación de rutas desde un teléfono celular.

La importación dependerá del formato de salida de la aplicación que estemos usando. Las aplicaciones de GPS de los celulares generalmente permiten exportar las rutas en formato .kml (GoogleMaps) o GPX. En general importar .kml es R es complicado, porque hay encontrar el identificador de la ruta. Importar .gpx es más simple y por lo tanto el método preferido. Para esto último usamos la función readGPX() del paquete plotKML. El formato .gpx es una especificación de xml para guardar rutas. De un archivo .gpx tiene cinco columnas: longitud, latitud, elevación, tiempo (año, mes, día, hora, minuto, segundo) y extensiones, que en la lista que regresa readGPX() están en $tracks.

El objetivo es generar un mapa raster de la zona CDMX en la que está registrada una ruta y etiquetar algunos de los puntos con el tiempo acumulado desde la partida a la llegada. Lo primero es muy fácil, lo segundo no tanto, pues el manejo de horas y fechas no es fácil en R.

library(plotKML)

## plotKML version 0.5-6 (2016-05-02)

## URL: http://plotkml.r-forge.r-project.org/

# Gráfico con inmutabilidad de datos. 
#====================================
#1. Importo, 
#2. formateo el timestamp, 
#3. Genero la suma acumulada de diferencia en minutos. 
#4. Selecciono una fila cada 50
#Grafico con etiquetas. 

#Posibles mejoras de código: 
#estudiar mutate lag() y encontrar una solución más elegante al problema de la suma acumulada. Básicamente el problema del NA al pcipio
#Un mejor sistema para seleccionar una fila cada 50 (o arbitrario)
#Una importación de mapas mejor que la de google maps.
#Encontrar un algoritmo para identificar paradas en semáforos. 


readGPX("19_feb._2017_10_54_56.gpx")$tracks[[1]][[1]][] %>% #Porque son 3 listas anidadas.
  mutate(
    time=as.POSIXct(                            #Porque mutate no acepta POSIXlt
      strptime(time,                            #Regresa POSIXlt
               format = "%Y-%m-%dT%H:%M:%OS")   #timestamp de gpx
               )
         ) %>% 
  mutate(acum=round(         
           cumsum(c(0,       #diff regresa length-1 y mutate() no acepta largos diferentes.  
             diff(time)      #Diferencia en x y x[-1]
                   )
                 )/60)) %>%  #Paso directo a minutos. 
           mutate(etiquetas=
                    ifelse(row_number(acum) %in%  seq(1,nrow(.), 50), #Selecciono 1 cada 30 filas. 
                           acum, 
                           NA)                                          #El resto a NA
                  ) -> ruta 

#Cargamos el mapa desde GoogleMaps. 
ggmap(get_map (location=c(mean(ruta$lon), mean(ruta$lat)), source="google", maptype="roadmap", zoom=12)) +
  geom_point(data=ruta, aes(x=lon, y=lat), size=.5) +  
  geom_label_repel(data=ruta, aes(x=lon, y=lat, label=etiquetas))

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=19.425538,-99.15993&zoom=12&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false

#Mala opción, googlemaps no permite especificar el rectángulo, sólo un un centro y el nivel de zoom. stamen y osm sí, pero el servidor de osm funciona mal y stamen está roto para ggmaps (regresa un .png y get_map espera un .jpg).

4 Funciones mágicas.

El paquete mxmaps, de Diego Valle, permite hacer mapas coropléticos de México con una sintaxis simplificada. Si la cartografía que queremos producir se apoya en datos incluidos en las bases que instala el propio paquete no se requiere casi ninguna especificación, simplemente señalar la columna de interés y luego pasar un segundo comando para generar el mapa. Es perfectamente posible extender la funcionalidad de mxmaps:: agregando datos a las bases estandar. Para eso es necesario conocer la estructura de las bases de datos incluidas y agregar las columnas adecuadas. Con str(df_mxstate) y str(mx_municipio) podemos explorar la estructura de la base Estatal y Municipal, las dos unidades territoriales admitidas por este paquete. df_mxstate tiene 32 filas de largo, una por cada Entidad Federativa, e incluye la columna region con la clave INEGI de entidad, del 01 al 32. df_mxmunicipio tiene 2457 filas, una por cada municipio, y en la columna region la clave INEGI de 5 dígitos para municipios. Usando estas claves podemos hacer un join para agregar datos adicionales desde otras bases con igual unidad de análisis. La sintaxis básica para producir mapas a nivel municipios es df_mxstate#value <- df_mxstate$variable_de_interés, donde variable de interés es la que controlará el parámetro de color. Señalamos la variable de interés asignándola en la columna $value. Luego llamamos a la función mxmunincipio_choroppleth(df_mxmunicipio), es decir, señalamos el data.frame incluido en la librería en el que previamente señalamos la variable de interés en $values. Para mapas a nivel Estados usamos el data.frame df_mxstate y mxstate_choropleth(). En ninguno de los casos es necesario especificar los polígonos, están incluidos por defecto y esa es la ventaja de mxmaps.

#Instrucciones para instalar el paquete, es necesario ejecutarlo sólo una vez. 
#Mucha atención a los errores que regresa el instalador. El paquete tiene muchas dependencias no declaradas y hay que buscar cual produce el error, instalar ese paquete y así que hasta que funciona. 
#advertencia: el paquete es alpha, es decir, muy temprano en la etapa de desarrollo y aparentemente no tiene desarrollo activo.  ¿abandonado?
if (!require(devtools)) {
    install.packages("devtools")
}
devtools::install_github('diegovalle/mxmaps')

library(mxmaps)
#Coropléticos a nivel estatal. 

df_mxstate$value <- df_mxstate$pop   #Especifico en $value la columna que controlará el color. 
mxstate_choropleth(df_mxstate)       #Produzco el gráfico. 

#Presentación alternativa: mapa hexagonal. 

mxhexbin_choropleth(df_mxstate)      #Los presenta como hexágonos uniformes. 

#A nivel municipal. 

df_mxmunicipio$value <- df_mxmunicipio$pop_male - df_mxmunicipio$pop_female  #Diferencia entre población masculina y femenina. 
mxmunicipio_choropleth(df_mxmunicipio)

#Extender la funcionalidad agregando datos no incluidos en mxmaps.
#=================================================================


#Preparar los datos para unirlos. 

marginacion$region <- marginacion$`<U+FEFF>CVE_MUN` #Retomo el data.frame marginacion, que ya está en mi entorno de trabajo. 
                                                    #Le agrego una columna llamada region (sin tilde) con la clave municipal.
                                                    #Así coincide con el nombre de la columma clave de mx_dfmunicipios
#Unir los datos. 

df_mxmunicipio <- 
  inner_join(df_mxmunicipio,                        #El documento maestro es df_mxmunicio
             filter(marginacion, AÑO=="2015"),      #Al que agrego los datos de 2015 de de la base del IM (los filtro directamente aquí)
             by="region")                           #region es la columna de empate, presente en ambos df.

#Verificar el resultado

df_mxmunicipio %>% str()    #Se agregaron correctamente las variables del IM. 

## 'data.frame':    2457 obs. of  41 variables:
##  $ region             : chr  "01001" "01002" "01003" "01004" ...
##  $ state_code         : chr  "01" "01" "01" "01" ...
##  $ state_name         : chr  "Aguascalientes" "Aguascalientes" "Aguascalientes" "Aguascalientes" ...
##  $ state_name_official: chr  "Aguascalientes" "Aguascalientes" "Aguascalientes" "Aguascalientes" ...
##  $ state_abbr         : chr  "AGS" "AGS" "AGS" "AGS" ...
##  $ state_abbr_official: chr  "Ags." "Ags." "Ags." "Ags." ...
##  $ municipio_code     : chr  "001" "002" "003" "004" ...
##  $ municipio_name     : chr  "Aguascalientes" "Asientos" "Calvillo" "Cosío" ...
##  $ pop                : int  877190 46464 56048 15577 120405 46473 53866 8896 20926 20245 ...
##  $ pop_male           : int  425731 22745 27298 7552 60135 22490 26693 4276 10197 9982 ...
##  $ pop_female         : int  451459 23719 28750 8025 60270 23983 27173 4620 10729 10263 ...
##  $ afromexican        : num  532 3 10 0 32 3 13 13 4 0 ...
##  $ part_afromexican   : num  2791 130 167 67 219 ...
##  $ indigenous         : num  104125 1691 7358 2213 8679 ...
##  $ part_indigenous    : num  14209 92 2223 191 649 ...
##  $ metro_area         : chr  "Aguascalientes" NA NA NA ...
##  $ value              : int  -25728 -974 -1452 -473 -135 -1493 -480 -344 -532 -281 ...
##  $ <U+FEFF>CVE_MUN    : chr  "01001" "01002" "01003" "01004" ...
##  $ CVE_ENT            : chr  "01" "01" "01" "01" ...
##  $ ENT                : chr  "Aguascalientes" "Aguascalientes" "Aguascalientes" "Aguascalientes" ...
##  $ CVEE_MUN           : chr  "001" "002" "003" "004" ...
##  $ MUN                : chr  "Aguascalientes" "Asientos" "Calvillo" "Cosío" ...
##  $ POB_TOT            : int  877190 46464 56048 15577 120405 46473 53866 8896 20926 20245 ...
##  $ VP                 : chr  "-" "-" "-" "-" ...
##  $ ANALF              : num  2.06 4.47 4.8 4.35 3.26 3.41 3.53 3.2 4.23 4.92 ...
##  $ SPRIM              : num  9.54 20.89 24.18 16.55 13.73 ...
##  $ OVSDE              : num  0.31 3.84 0.55 1.24 0.44 1 1.97 3.51 1.85 4.89 ...
##  $ OVSEE              : num  0.16 1.04 0.41 0.57 0.37 0.42 0.52 1.32 0.63 1.97 ...
##  $ OVSAE              : num  0.72 1.17 0.86 0.47 0.73 0.9 1.52 1.23 0.94 2.35 ...
##  $ VHAC               : num  18 30.6 30.5 34.6 27 ...
##  $ OVPT               : num  0.63 1.69 0.93 1.42 1.03 0.62 0.91 0.44 1 1.65 ...
##  $ PL<5000            : num  8.73 100 50.76 100 45.17 ...
##  $ PO2SM              : num  31.1 52.8 62 49.8 33.8 ...
##  $ OVSD               : chr  "-" "-" "-" "-" ...
##  $ OVSDSE             : chr  "-" "-" "-" "-" ...
##  $ IM                 : num  -1.676 -0.565 -0.698 -0.674 -1.256 ...
##  $ GM                 : chr  "Muy bajo" "Bajo" "Bajo" "Bajo" ...
##  $ IND0A100           : chr  "-" "-" "-" "-" ...
##  $ LUGAR_NAC          : chr  "2 408" "1 669" "1 799" "1 773" ...
##  $ LUGAR_EST          : int  11 1 5 4 10 8 7 6 2 3 ...
##  $ AÑO                : chr  "2015" "2015" "2015" "2015" ...

df_mxmunicipio$value <- df_mxmunicipio$GM     #Ya estoy usando una variable importada del segundo data.frame

mxmunicipio_choropleth(df_mxmunicipio) +      #Llamada usual de mxmunicipio_choropleth
  labs(title="Grado de marginación municipal", #Agrego sintaxis de ggplot2 para etiquetado
       subtitle="Mapeado con mxmaps",
       caption="Datos estadísticos CONAPO\nDatos geográficos ¿?\nProyección azimutal de áreas equivalentes.") +
  coord_map("azequalarea") +                 #Proyección del mapa
  scale_fill_manual(values=
                      c("red", "orange", "yellow", "red4", "green")) #Escala discreta de colores.

## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

5 Librerías requeridas.

library(rgdal)     #Para importar shapefiles. 
library(knitr)     #Para formatear tablas
library(mapproj)   #Para cambiar las proyecciones. 
library(ggrepel)   #Para controlar etiquetado.
library(ggmap)     #Para importar rasters y producir mapas. 
library(plotKML)   #Levanta muchísimas dependencias en la instalación.
library(mxmaps)    #devtools::install_github('diegovalle/mxmaps')
library(tidyverse)

---

1. mpaladino@mora.edu.mx↩